Unitary circuits for strongly correlated fermions 
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We introduce a scheme for efficiently describing pure states of strongly correlated fermions in higher di- 
mensions using unitary circuits featuring a causal cone. A local way of computing local expectation values 
is presented. We formulate a dynamical reordering scheme, corresponding to time-adaptive Jordan-Wigner 
transformation, that avoids nonlocal string operators. Primitives of such a reordering scheme are highlighted. 
Fermionic unitary circuits can be contracted with the same complexity as in the spin case. The scheme gives 
rise to a variational description of fermionic models not suffering from a sign problem. We present numerical 
examples on 9 x 9 and 6x6 fermionic lattice model to show the functioning of the approach. 



Strongly correlated quantum lattice models still pose some 
of the most intriguing problems in quantum physics, presum- 
ably being at the basis of phenomena such as high-temperature 
superconductivity [l]. To describe such models theoretically 
is often very difficult, an obvious obstacle being the pro- 
hibitive dimension of Hilbert space even for moderately sized 
systems when naively representing ground states. In recent 
years, it has increasingly become clear, however, that typi- 
cal ground states of physically meaningful local models are 
lurking in some corner of Hilbert space, one that can often 
even be identified ('Z'-H^. Hence to faithfully describe such a 
state, even though it may be impossible to parametrize the en- 
tire physical space, one merely has to parametrize those quan- 
tum states from the relevant set, requiring significantly fewer 
parameters. The most prominent and, to date, most power- 
ful incarnation of this idea is provided by the density-matrix 
renormalization group (DMRG) method l2l[3l, requiring only 
linearly many parameters in the system size, but still giving 
an exceptionally good description of gapped one-dimensional 
spin chains and a very reasonable one for critical chains. 

To reaUze such an idea in higher-dimensional systems is 
significantly more difficult, although progress has been made 
in recent years ||31 4T3]| . specifically when it comes to generaliz- 
ing DMRG ideas to higher dimensions, including variations of 
tensor product states or projected entangled pair states, vari- 
ants with graph enhancement, or the promising multiscale en- 
tanglement renormalization (MERA). This is a scale-invariant 
approach related to renormalization that has been proposed in 
Ref. |6| and implemented in Refs. irT UTTl . Area laws 1 12| of- 
ten point to the part of Hilbert space that is being occupied. 

To describe higher dimensional fermionic systems, such as 
the Fermi-Hubbard model itself [l], in such a fashion, ap- 
pears particularly promising, but also particularly challenging. 
Here, other key methods of describing quantum lattice models 
like the powerful quantum Monte Carlo method are hampered 
by the sign problem 1 14|. Surely, fermionic models can read- 
ily be represented as spin models, yet at the expense of losing 
locality (or by increasing the locality region of Hamiltonians 
ifTSl ). If one considers the term fjfk, j < k, then its spin rep- 



resentation under the Jordan-Wigner transformation (JWT) is 
(Jj ® at ® cr+. 

j<l<k 

For a pair of nearest-neighbor sites {j, k) on a d x d lattice, 
the occurring string operator that is supported not only on the 
spins associated with j and k, but in fact on all spins between 
j and fc, will typically have a length linear in d. It should be 
clear that no order can be chosen to let this apparent problem 
disappear. 

In this work, we present a method for studying strongly cor- 
related fermionic models using quantum circuits, that is, cir- 
cuits of fermionic gates, in a way that is not overburdened by 
string operators: When describing the system and computing 
local expectation values, one has to deal with operators hav- 
ing support identical to that of a coiTesponding spin system 
(fermionic gates of the circuit replaced by regular ones), and 
strings can be made to disappear. The key point is to acknowl- 
edge that while any fixed order will give rise to the aforemen- 
tioned problem, we are not necessarily forced to pick any or- 
der in the first place. Instead, one can employ a dynamical 
reordering of the fermionic modes in the essential part of the 
lattice, reordering and projecting out particles "on the fly," in 
dynamical JWTs. We show that this can be consistently done, 
not altering expectation values of parity-preserving operators. 
In this way, we find that to describe fermionic lattice models 
is, in this sense, just as hard or easy as describing spin models. 
More formally, the contraction complexity of the circuit is the 
same. 

Fermions. We first prepare the ground of the dynamical re- 
ordering idea. The algebra GiL) of a set of n fermionic modes 
is spanned by products of anticommuting fermionic operators 
: J = 1, . . . , n}, with M = 0, {f,jl} - S,,k. 
Physical operators form the subalgebra J-{L), the so-called 
physical algebra, of operators respecting the fermion number 
parity |16, 17|. In practical terms, this means that they are 
even polynomials in the fermionic creation and annihilation 
operators. Hence, the physical algebra splits into a direct sum 
of an even and an odd part, F{L) = F^°'^'^\L). 

For a subset / C L of some sites, one similarly finds the phys- 
ical algebra J-{I). 
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Fermionic unitary circuits and multiscale entanglement 
renormalization. A fermionic unitary is a parity-preserving 
unitary gate acting on fermionic modes, U = exp{iH), where 
H S J^{L) is a hermitian operator. In a circuit |6 -9|, unitaries 
will typically not have support on the entire lattice, but will be 
local; that is, it will have a small support / independent of 
n such that U E ^{I)- A fermionic unitary circuit is an or- 
dered product of fermionic unitaries. For most of this work, 
we in fact consider conjugation, in particular for the evalu- 
ation of expectation values of local observables A e J^iL), 
with respect to states that have been prepared by applying a 
fermionic unitary circuit to the vacuum 1 0) , 

E = {0\...U2UiAUlu^^...\0), (1) 

or by applying it to some other with odd fermion number par- 
ity, for example, |0) instead of |0) . In fact, there is a natural 
discrete time label t in a circuit, in that t = 0, . . . ,T labels 
the time at which a single gate is being applied. What is more, 
we focus on circuits that feature a causal cone (see Fig. [TJ, 
most prominently present in the MERA in the original sense 
of Ref. |6|: The cone is formed by those unitaries in a cir- 
cuit that cannot be sequentially canceled with their conjugates 
from the dual vector in Eq. ([T]i due to the presence of a local 
observable A that is supported only on a few sites, typically 
nearest neighbors. Efficient schemes will have a causal cone 
of a fixed width, in that a local operator wiU only touch a con- 
stant number of unitaries for each time step. The sequence of 
sequentially computing the expectation values respecting the 
discrete time order gives rise to a contraction. We can also al- 
low for tensors different from unitary gates and partial traces. 
We also allow time steps where partial projections are applied 
onto the fermionic vacuum in some mode. 

We now turn to aspects that are specific for fermions. A 
first key property from the elements of disjoint subalgebras 
U e J"(/) and V G J^{J), I Ci J = ^, respecting parity, is 
that their elements commute, [U, V] — 0. This has the im- 
portant consequence that just as for spin systems, all unitaries 
outside a causal cone can be canceled, and the fermionic vari- 
ant inherits the same cone as the spin system. We now see that 
the computation of the expectation value Eq. ([T]) can be done 
without having to deal with string operators outside the cone. 

Jordan-Wigner transformations. Fermionic operators are 
encoded in an occupation number representation, a repre- 
sentation necessarily depending on the chosen order of the 
fermionic modes. As we will later make updates with dif- 
ferent fermionic orderings at different times, it makes sense 
to first highlight the status of the JWT. For the entire sys- 
tem L one can take some order O corresponding to an el- 
ement of the symmetric group Sn, for example the trivial 
order O — (1, . . . , n). For this order O, a basis for the 
Hilbert space is given in the occupation number representa- 
tion, . . . ,i„)o = ifoiY' ■ ■ ■ (/o„)'"|0)' giving an iden- 
tification of the total Hilbert space with the (C^)*^", the 
Hilbert space of n spins 1/2. Expressing fermionic operators 
in the occupation number representation amounts to a map 




Figure 1 : Example for the evaluation of a local observable (black 
box) with respect to a fermionic MERA (gray boxes). Each gray box 
represents a single fermionic unitary. Dark gray boxes make up the 
causal cone of the observable; all others lie outside of it and cancel 
each other. Only for conceptual clarity, a ID MERA is depicted, not 
a 2D one, for which dynamical reordering becomes relevant. 

Jo ■ G{L) M{L), where M{L) is the algebra of linear 
operators on ((D^)^". This map is the well-known JWT, de- 
pending on L and the chosen order O € Sn of the fermionic 
modes; for example, 

^o(/oj = (0^r) ®a+ (2) 

where cr^'^'^ denote the familiar Pauli operators and cr^ — 
(crj + ia^)/2. In this light, a JWT simply gives fermionic 
operators in the occupation number (or spin) representation 
for a given ordering O. We hence identify local fermionic 
operators with non-local operators in the respective spin sys- 
tems. Also, physical operators contain strings that depend on 
the order chosen. Consider a local physical operator acting 
on sites /, the spin representation will in general also be sup- 
ported on the spins in between with respect to order O. For 
example, for the local fermionic number operator one finds 
Jo ifo foj ) ~ '^j' the occupation number representation 

of /o,/o. 

Mfljo,)=<Jj® (g) 'J!®a+ (3) 

j<l<k 

with j < k, contains a string of Pauli operators on [j, k]. 

For operators A E J^{L) parity conservation manifests it- 
self in the occupation number representation for an O as fol- 
lows: Splitting the occupation number basis into two parts, 
depending on the parity of ik, the spin representation 
a = Jo{A) of A attains the form a = a<'=™"' © a^°^'^\ 

Dynamical reordering. We now describe a new scheme to 
contract a fermionic unitary circuit, henceforth called dynam- 
ical reordering, and demonstrate that the contraction can be 
done with essentially the same complexity compared to the 
corresponding spin circuit: We take the "time order" in con- 
tracting, with t — corresponding to the local operator A 
itself In turn, in our fermionic variant of the MERA |6|, the 
top layer with largest t defines the parity, similar to fixing the 
topological degrees of freedom in Kitaev's toric code |18|. 
The general idea of dynamical reordering is to order only the 
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modes we are currently using at a time i in a nontrivial fashion 
with order O* expHcitly dependent on t. As we proceed, new 
modes appear in the causal cone, and they are added to the 
description. When they are discarded due to projection, they 
are removed from the description O* as t evolves. Hence, we 
arrive at an entirely local description of the fermionic tensor 
network, using the subsequent rules (see Fig.[2]l. 

(a) Local representation of fermionic local operators: Con- 
sider a local operator A e ^(/); that is, it is an even poly- 
nomial in the operators {/j, /|} for j £ I. The spin repre- 
sentation of the fermionic operator A in order O is given by 
a = Jo{^)- This representation only involves k = \I\ spins. 

(b) Reordering fennionic modes: Consider some fermionic 
operator A and an order O* at time t. The new spin rep- 
resentation at a time t + 1 in the new order 0*+^ reads 
Jot+i{A) = pJotiA)p, with p = Hi 1 ® -^01,0* ^ ® 1' 
soiol^, = |0,0)(0,0| + |0,l)(l,0| + |l,0)(0,l|-il,l)(l,l|, 
where the i are chosen such that the corresponding sequence 
of nearest-neighbor pair permutations of modes gives rise to 
the overall permutation 0*+^ = 7r(0*); so familiar fermionic 
swaps can be applied to the local spin representation. 

(c) Prepending of fermionic modes: Let us have an operator 
A at some time t, represented as a = Jq* (A) in the order 
O*. We now just prepend a new fermionic mode, yielding the 
new order 0*+^ — {k, O*). Then the new representation of 
a' = Jot+^ (^) is, in the even and odd sectors, given by 

(^.<)(evenorodd) ^ |q^^q| ^ ^(even or odd) _^ |-^^^-^| ^ ^(odd or even) 

( d) Conjugation with fennionic unitaries having the same 
support: Consider a subset of modes / at time t. Let 
U,A£ J- {I) be fermionic operators with support in /. Then 
Jot+i{UAU'^) = Jot (U) Jot (A) Jot {W) if the same order 
0*+^ = O* is taken. If U is stored in a different order, the 
permutation rule has to be applied first. 

(e) Partial trace over the first mode: Let us have some op- 
erator A at time t in order O*, and we trace out the first mode 
0{; that is, O* = {0{, 0*+^). Then the spin representation of 
the resulting fermionic operator is, with k — |0*|, 

Jot+i(ti-otA) = ^{i,j2,---,jk\Jot{A)\i,i2,---,ik) 

X |j2, ■ • ■ , jfc)(«2, ■ • ■ , jfcl- (5) 

(f) Partial projection over the first mode: Let us consider 
the same orderings O* and 0*+^ as for rule (e), but instead 
of tracing out mode r ~ 0{ from operator A, we now want 
the projection of the mode to the empty state or that of a filled 
mode. The corresponding projectors in J-{I) are P,-*'^ = /,./^ 
and Pr^^ = flfr, respectively, so that the projected fermionic 
operators are Pr''^ APr^\ i = 0, 1. Applying this to Eq. (|5]l, 
we find 

Jot^i{X^rP^APi^) - 5] |J2, . . ■,3k){i2. ..-M 

X {i,32, ■ ■ ■ ,jk\Jot{A)\i,i2,. . ■ ,ik)- (6) 




Figure 2: The basic primitives of (a) representing local operators, (b) 
reordering fermionic modes, (c) adding fermionic modes, (d) con- 
jugating with unitaries, (e) doing partial traces, and (f) generating 
partial projections. 

The main point now is that if one follows this recipe of 
rules (a)-(f), and keeps only local descriptions at each time, 
involving the nontrivial support, one gets a value identical 
to Eq. ([T]i as if one had performed a global JWT and com- 
puted it-then with a significant overhead in complexity. To 
prove the reordering rule for two modes with labels j, k with 
j = Ol and k — 0,*^i, that is, two modes neighbor- 
ing in the order O*, we need to show Jo* (S'j j.) = 
Sj\fc-^o'(^)s].fe = Jo*+^{A), where Sj^k acts as Sj^kfjS],^ = 
fk and Sj^kfkS]j^, = fj. One finds S,.k - 1 - fjfj - 

flfk + f]fk + flfj MM- The reordering rule for more 
than two modes follows by iteration. To derive the rule of 
adding modes, note that the local JWT gives Jotifoj,) — 
((XjfjTi^crf) «) cT^, and hence we obtain for (9*+^ = {x,0*) 
that Jot+i{fot+i) = Jo'+i(/ot) = erf «) Jo'(/of)- Ex- 

fc + 1 

ploiting the conservation of the fermion number parity for 

any A e Jo*+i(^) = E.K)'" ® Jo'(A.) = 

1 (8i Jot (A), where collects all monomial terms of de- 
gree 2iy such that A = With a' = Jot+i{A) 
and a = Jot (A), for the even and odd sectors, this 
gives Eqs. Q. To show Eq. (|5]l, note that for any op- 
erator B not acting on mode Oi, that is, B E F{I \ 
{Ol}), matrix elements obey (ii, . . . , ik\oB\ii, . . ■,jk)o = 
(«2, ■ ■ ■,ik\o'B\j2, ■ ■ ■,jk}o', where O = (Oi,0') and 
k ^ \0\. With 

troi(^) = '\3\^31 ■ ■ ■ ,jk\oA\ji,i2, ■ ■ ■,ik)o 

X Ih, ■ ■ ■,3k)o'{i2, ■ ■ -Jklo', (7) 

it follows that tr{AB) = tr[troi(^)S] for any A e J^{I)- 
This means that troi is the fermionic partial trace for 
mode Ol. The resulting occupation number representation 
Jo' {iTOi A) of the partial trace is given in Eq. (|5]). Rules (c), 
(e), and (f) were given for the case where only the very first 
mode Ol of the ordering is affected. Combining those rules 
with the rule for reordering modes, they are generalized to the 
case where an arbitrary mode Oj is affected. In the corre- 
sponding expressions this results in extra sign factors. For the 
generalization of rule (d), note that when one wants to perform 
a contraction in which the support is not the same for the two 
involved operators, one can add the corresponding modes in 
both operators, reorder and apply (d). With theses basic build- 
ing blocks, we can thus contract fermionic tensor networks 
and also build, for example, reduced density matrices fSJ. 
Observation (Computing expectation values). When com- 
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Figure 3: Numerical examples of the ground-state energy _Eo for 
several 9x9 and 6x6 (boxes and triangles) instances of strongly cor- 
related fermionic models, H = iY,{j.k)U] fk + h.c.) + fjfj + 

fe> fjfifkf''' f*^^^ models u = are marked. As the 
MERA ansatz is a truly variational ansatz, providing upper bounds 
to the energy, success is guaranteed by giving lower bounds as pro- 
vided by Anderson bounds (lines). This is based on an exact diag- 
onalization of 5 x 5 fermionic models, taking care of arising string 
operators. Even for this lowest order 2D MERA, in which, in each 
renormalization step, nine fermionic modes are mapped to one, we 
find an impressive precision. The inset shows the convergence in 
time for a 9 X 9 free model (u = 0, t = 0.3) in a minimization using 
imaginary time evolution. In the same way, correlators and on-site 
properties can be efficiently computed. 



puling expectation values E = {<I>\ . . . U2U1HUIUI . . . |*) 
where U = U1U2 ■ ■ ■ is a fermionic circuit on an elemen- 
tary state vector |\['), then the support of the causal cone at 
each time step is exactly the same as if one had a spin sys- 
tem with the same topology. Hence, a fermionic circuit acting 
on a higher-dimensional fermionic system can be described 
entirely locally with the same time and memory complexity. 

This observation renders known algorithms useful even in 
the case of strongly correlated fermions, with little modifi- 
cation. Notably, for MERA, the computational effort for the 
computation of a local expectation value scales as 0(log I) in 
an cubic lattice and polynomially in the refinement param- 
eter |6 |, hence giving rise to an indeed efficient algorithm, up 
to, at most, a small constant the same as for bosons or spins. 
It is noteworthy that the entire machinery developed for spin 
systems can, with slight modification, be used here. The pre- 
ceeding observation is independent of whether the system is 
free or strongly correlated. As it turns out, general fermionic 
gates also can be taken account, by directly contracting appro- 
priately ordered fermionic operators. One can show that the 
computational requirements remain unchanged giving rise to 
a general theory of fermionic tensor networks. 

Time evolution. One also finds the following: For the time 
evolution, both real and imaginary, of local fermionic Hamil- 
tonians, we aim at determining updates of each the fermionic 
unitaries forming the circuit such that )^ (?72)^ . . . is an 
optimal approximation to the evolved exp{xH)ulU2 ■ ■ ■ ji/'), 
again with the same complexity as for a qubit system. 

Numerical implementation. We have realized a scal- 
able numerical implementation of this idea, applied to two- 



dimensional fermionic lattice models. To show the function- 
ing of this approach as a proof of principle, we present bench- 
mark examples in Fig.|3] 

Summary. In this work, we have introduced a notion of 
time-adaptive JWTs allowing us to contract fermionic uni- 
tary circuits with the same complexity as for the correspond- 
ing spin model. This opens up a way to efficiently describe 
fermionic higher-dimensional models without intrinsic algo- 
rithmic problems, such as the sign problem in quantum Monte 
Carlo sampling. It is the hope that this work stimulates further 
interest in numerically assessing strongly correlated fermionic 
models using unitary circuits. 

Note added: We acknowledge discussions with U. Scholl- 
woeck. We also warmly acknowledge an early key discus- 
sion with F. Verstraete, and both this work and Ref. [20 1 were 
sparked off from this initial discussion. In the meantime, fur- 
ther approaches to the problem have appeared ||2T| . 
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